#################
#Data management#
#################
#orginal dataset "Marokko 2008.csv" --> one row per department, no department with same identifier
#merging of departments so that equal size
library(psych)
mydata<-read.table("Data.csv", header=TRUE, sep=",")
describe(mydata)

table(mydata$ministry,mydata$department) 

mydata$department<- ifelse(mydata$ministry ==4 & mydata$department == 1, 0, mydata$department)
mydata$department<- ifelse(mydata$ministry ==4 & mydata$department == 6, 7, mydata$department)
mydata$department<- ifelse(mydata$ministry ==6 & mydata$department == 7, 0, mydata$department)
mydata$department<- ifelse(mydata$ministry ==6 & mydata$department == 8, 9, mydata$department)
mydata$department<- ifelse(mydata$ministry ==7 & mydata$department == 22, 0, mydata$department)
mydata$department<- ifelse(mydata$ministry ==7 & mydata$department == 10, 9, mydata$department)
mydata$department<- ifelse(mydata$ministry ==7 & mydata$department == 12, 0, mydata$department)
mydata$department<- ifelse(mydata$ministry ==7 & mydata$department == 11, 10, mydata$department)
mydata$department<- ifelse(mydata$ministry ==8 & mydata$department == 16, 0, mydata$department)
mydata$department<- ifelse(mydata$ministry ==8 & mydata$department == 14, 7, mydata$department)
mydata$department<- ifelse(mydata$ministry ==8 & mydata$department == 21, 0, mydata$department)
mydata$department<- ifelse(mydata$ministry ==8 & mydata$department == 24, 8, mydata$department)

######################
#Non-imputed data set#
######################

##Operationalization Dependent Variable##
library(psych)
mydata$DemGov1<-(mydata$item1 + mydata$item2 + mydata$item3 + mydata$item4 + mydata$item5 + mydata$item6 + mydata$item7 + mydata$item8)
mydata$DemGov2<-(mydata$item1 + mydata$item2 + mydata$item3 + mydata$item5 + mydata$item6 + mydata$item7 + mydata$item8)
mydata$DemGov3<-(mydata$itemN1 + mydata$itemN2 + mydata$itemN3)
mydata$DemGov4<-(mydata$DemGov1 + mydata$DemGov3)
psych::describe(mydata$DemGov1)
psych::describe(mydata$DemGov2)
psych::describe(mydata$DemGov3)
psych::describe(mydata$DemGov4)

#Table 2. Summary of items statements
psych::describe(mydata$item1)
psych::describe(mydata$item2)
psych::describe(mydata$item3)
psych::describe(mydata$item4)
psych::describe(mydata$item5)
psych::describe(mydata$item6)
psych::describe(mydata$item7)
psych::describe(mydata$item8)
psych::describe(mydata$itemN1)
psych::describe(mydata$itemN2)
psych::describe(mydata$itemN3)

#reliability of scale#
names(mydata)
DemGov1<-mydata[,c(25,26,27,28,29,30,31,32)]
DemGov2<-mydata[,c(25,26,27,29,30,31,32)]
DemGov3<-mydata[,c(34,35,33)]
DemGov4<-mydata[,c(25,26,27,28,29,30,31,32,34,35,33)]
psych::alpha(DemGov1)
psych::alpha(DemGov2)
psych::alpha(DemGov3)
psych::alpha(DemGov4)

##Operationalization Treatment Variables##
#International education: frequency#
table(mydata$abroad_EU)
table(mydata$abroad_USA)
table(mydata$abroad_USA & mydata$abroad_EU)

#International education: differences in attitude betw. categories
#levels=c(0,1,2,3)
#labels=c("No stay abroad", "stay abroad only in the US", "only in Europe", "in both Europe and USA")

abroad_NO <- ifelse (mydata$abroad_EU==0 & mydata$abroad_USA==0,0,0) 
abroad_US <- ifelse (mydata$abroad_EU==0 & mydata$abroad_USA==1,1,0)
abroad_EUR <- ifelse (mydata$abroad_EU==1 & mydata$abroad_USA==0,2,0)
abroad_EUUS <- ifelse (mydata$abroad_EU==1 & mydata$abroad_USA==1,3,0)
 
group <- (abroad_NO+abroad_US+abroad_EUR+abroad_EUUS)
table(group)
value1<-(mydata$item1+mydata$item2+mydata$item3+mydata$item4+mydata$item5+mydata$item6+mydata$item7+mydata$item8)/8
value2<-(mydata$item1+mydata$item2+mydata$item3+mydata$item4+mydata$item5+mydata$item6+mydata$item7+mydata$item8+mydata$itemN1+mydata$itemN2+mydata$itemN3)/11
mydata2<-data.frame(group, value1, value2)

bartlett.test(value1~group, data=mydata2) 
bartlett.test(value2~group, data=mydata2) 
#p > 0.05; we can assume that the variances are equal

aov.abr1 = aov(value1~group, data=mydata2)
aov.abr2 = aov(value2~group, data=mydata2)
summary(aov.abr1)  
summary(aov.abr2)                                   
print(model.tables(aov.abr1,"means"),digits=3)
print(model.tables(aov.abr2,"means"),digits=3)
boxplot(value1~group)
boxplot(value2~group)

#Networks: frequency#
mytable.Twin<-table(mydata$netw_Twin)
mytable.EUMS<-table(mydata$netw_EUMS)
mytable.UNDP <-table(mydata$netw_UNDP)
mytable.USA <-table(mydata$netw_USA)
mytable.WB <-table(mydata$netw_WB)
mytable.Japon <-table(mydata$netw_Japon)
prop.table(mytable.Twin) * 100
prop.table(mytable.EUMS) * 100
prop.table(mytable.UNDP) * 100
prop.table(mytable.USA) * 100
prop.table(mytable.WB) * 100
prop.table(mytable.Japon) * 100

#netw = Twinning program, Alt = alle anderen Programme)
mydata$netw<-ifelse(mydata$netw_Twin==1,1,0)
mydata$netw <-ifelse(mydata$ministry==11,0,mydata$netw)

mydata$Alt<-ifelse(mydata$netw_Twin==0 & (mydata$netw_EUMS==1 | mydata$netw_UNDP==1 |
mydata$netw_USA==1 | mydata$netw_WB==1 | mydata$netw_Japon==1), 1, 0)

mytable.netw<-table(mydata$netw)
prop.table(mytable.netw) * 100

mytable.Alt<-table(mydata$Alt)
prop.table(mytable.Alt) * 100

mytableA <- table(mydata$ministry,mydata$netw) 
mytableB<-table(mydata$ministry, mydata$Alt)
mytableA
mytableB
 
#creation of alternative treatment & control variables
mydata$abroad <- ifelse(mydata$abroad_EU == 1 | mydata$abroad_USA == 1,1,0)
mydata$infl<-ifelse(mydata$abroad==1|mydata$netw==1|mydata$media==1,1,0)
mydata$polit <- ifelse(mydata$ministry==2 | mydata$ministry==7 | mydata$ministry==8 |
mydata$ministry==9 | mydata$ministry==10 | mydata$ministry==11,1,0)

#intercorrelation internat. education and foreign media#
mytable<-table(mydata$abroad, mydata$media)
mytable
chisq.test(mytable)


###Descriptive Statistics###
#Table Appendix: Basic Information of Variables
library(Hmisc)
data.link<-data.frame(mydata$netw, mydata$media, mydata$abroad)
names(data.link)<-c("netw", "media", "abroad")
rcorr(as.matrix(data.link), type="spearman")

library(psych)
describe(mydata$DemGov4)
describe.by(mydata$DemGov4, mydata$infl)
describe.by(mydata$DemGov4, mydata$netw)
describe.by(mydata$DemGov4, mydata$media)
describe.by(mydata$DemGov4, mydata$abroad)

#gender
describe(mydata)
table(mydata$genderF)

library(prettyR)
describe(mydata$netw)
freq(mydata$netw)
describe(mydata$abroad)
freq(mydata$abroad)
describe(mydata$media)
freq(mydata$media)


describe(mydata$Alt)
freq(mydata$Alt)
describe(mydata$polit)
freq(mydata$polit)
describe(mydata$age)
freq(mydata$age)
describe(mydata$genderF)
freq(mydata$genderF)
describe(mydata$position)
freq(mydata$position)
describe(mydata$educ_level)
freq(mydata$educ_level)
describe(mydata$subject1)
freq(mydata$subject1)
describe(mydata$subject2)
freq(mydata$subject2)
describe(mydata$subject3)
freq(mydata$subject3)
describe(mydata$lang_e)
freq(mydata$lang_e)
describe(mydata$lang_f)
freq(mydata$lang_f)
describe(mydata$educ_subject)
freq(mydata$educ_subject)

#multicollinearity#
#variance inflation factor (VIF) quantifies the severity of multicollinearity: if values are less than 10, no cause for concern; average of VIF values no greater than 1, no cause for concern#
library(car)
modelVIF1<-lm(mydata$DemGov4~mydata$media+mydata$abroad+mydata$netw+mydata$Alt+mydata$polit+mydata$position+mydata$educ_level+mydata$subject1+mydata$subject3+mydata$lang_e+mydata$genderF+mydata$age, na.action=na.exclude)
modelVIF2<-lm(mydata$DemGov4~mydata$media+mydata$abroad+mydata$netw, na.action=na.exclude)
vif(modelVIF1)
vif(modelVIF2)

#Figure 1: Beanplots##
library(lattice)
library(beanplot)

mydata$Infl <-with(mydata, ifelse (infl==1,"A 1", "A 0"))
mydata$Netw <-with(mydata, ifelse (netw==1,"D 1", "D 0"))
mydata$Media <-with(mydata, ifelse (media==1, "C 1", "C 0"))
mydata$Education <- with(mydata, ifelse (abroad==1, "B 1", "B 0"))

data.infl <- mydata[,c("DemGov4", "Infl")]
data.netw <- mydata[,c("DemGov4", "Netw")]
data.media <- mydata[,c("DemGov4", "Media")]
data.abroad<- mydata[,c("DemGov4", "Education")]
colnames(data.infl)[2] <- "Types"
colnames(data.netw)[2] <- "Types"
colnames(data.media)[2] <- "Types"
colnames(data.abroad)[2] <- "Types"
data.beans <- rbind(data.infl, data.netw, data.media, data.abroad)

pdf("figure1.pdf", height=8)
par (mai = c(0.6, 0.7, 0.7, 0.4))
beanplot(DemGov4 ~ Types, data = data.beans, side = "both",ll=0.0,
what=c(1,1,1,1),log="",bw=0.75,
border = NA, col = list("black", c("grey", "white")),
names=c("Influences", "Education", "Media", "Networks"),
cut=0.5, beanlinewd=0.02)
legend("bottomleft", fill = c("black", "grey"),
legend = c("Comparison group", "Treatment group"))
dev.off()

##Distribution of Sample##

#Table A1
#levels=c(1,2,3,4,6,7,8,9,10,11)
#labels=c("Health", "Economic and General Affaires", "Housing, Urban Development and Planning", 
#"National Education, Higher Education, and Scientific Research", "Energy, Mining, Environment, and Water", 
#"Economy and Exterior Finances", "Agriculture, Rural Development, and Fishing",
#"Equipment and Transport","Industry, Trade, and New Technologies", "Foreign Commerce")
mytable1 <- table(mydata$ministry,mydata$netw) 
mytable2<-table(mydata$ministry, mydata$abroad)
mytable3<-table(mydata$ministry, mydata$media)
mytable1
mytable2
mytable3

mytable1.n <- table(mydata$position,mydata$netw)
#levels=c(0,1)
#labels=c("Administrative Staff", "Directorate")
mytable2.n<-table(mydata$educ_level, mydata$netw)
#levels=c(0,1)
#labels=c("University diploma", "Postgraduate")
mytable3.n<-table(mydata$educ_subject, mydata$netw)
#levels=c(1,2,3)
#labels=c("Law / Economics", "Natural Sciences", "Public Administration / Political Science")
mytable4.n<-table(mydata$lang_e, mydata$netw)
mytable5.n<-table(mydata$lang_f, mydata$netw)
mytable6.n<-table(mydata$genderF, mydata$netw)
mytable1.n
mytable2.n
mytable3.n
mytable4.n
mytable5.n
mytable6.n


mytable1.a <- table(mydata$position,mydata$abroad) 
mytable2.a<-table(mydata$educ_level, mydata$abroad)
mytable3.a<-table(mydata$educ_subject, mydata$abroad)
mytable4.a<-table(mydata$lang_e, mydata$abroad)
mytable5.a<-table(mydata$lang_f, mydata$abroad)
mytable6.a<-table(mydata$genderF, mydata$abroad)
mytable1.a
mytable2.a
mytable3.a
mytable4.a
mytable5.a
mytable6.a

mytable1.m <- table(mydata$position,mydata$media) 
mytable2.m<-table(mydata$educ_level, mydata$media)
mytable3.m<-table(mydata$educ_subject, mydata$media)
mytable4.m<-table(mydata$lang_e, mydata$media)
mytable5.m<-table(mydata$lang_f, mydata$media)
mytable6.m<-table(mydata$genderF, mydata$media)
mytable1.m
mytable2.m
mytable3.m
mytable4.m
mytable5.m
mytable6.m

##Regression Analyses##
#Independent effects#
#subject1=law/econ; subject2 = natural sciences; subject3 = public admin

#positively-oriented items, i.e. DemGov1#
model1P<-lm(mydata$DemGov1~mydata$media+mydata$Alt+mydata$polit+mydata$position+mydata$educ_level+
mydata$subject1+mydata$subject3+mydata$lang_e+mydata$genderF+mydata$age, na.action=na.exclude)
model2P<-lm(mydata$DemGov1~mydata$abroad+mydata$Alt+mydata$polit+mydata$position+mydata$educ_level+
mydata$subject1+mydata$subject3+mydata$lang_e+mydata$genderF+mydata$age, na.action=na.exclude)
model3P<-lm(mydata$DemGov1~mydata$netw+mydata$Alt+mydata$polit+mydata$position+mydata$educ_level+
mydata$subject1+mydata$subject3+mydata$lang_e+mydata$genderF+mydata$age, na.action=na.exclude)
model4P<-lm(mydata$DemGov1~mydata$media+mydata$abroad+mydata$netw+mydata$Alt+mydata$polit+
mydata$position+mydata$educ_level+mydata$subject1+mydata$subject3+mydata$lang_e+mydata$genderF+mydata$age, na.action=na.exclude)
summary(model1P)
AIC(model1P)
summary(model2P)
AIC(model2P)
summary(model3P)
AIC(model3P)
summary(model4P)
AIC(model4P)

#negatively-oriented items, i.e. DemGov3#
model1N<-lm(mydata$DemGov3~mydata$media+mydata$Alt+mydata$polit+mydata$position+mydata$educ_level+
mydata$subject1+mydata$subject3+mydata$lang_e+mydata$genderF+mydata$age, na.action=na.exclude)
model2N<-lm(mydata$DemGov3~mydata$abroad+mydata$Alt+mydata$polit+mydata$position+mydata$educ_level+
mydata$subject1+mydata$subject3+mydata$lang_e+mydata$genderF+mydata$age, na.action=na.exclude)
model3N<-lm(mydata$DemGov3~mydata$netw+mydata$Alt+mydata$polit+mydata$position+mydata$educ_level+
mydata$subject1+bmydata$subject3+mydata$lang_e+mydata$genderF+mydata$age, na.action=na.exclude)
model4N<-lm(mydata$DemGov3~mydata$media+mydata$abroad+mydata$netw+mydata$Alt+mydata$polit+
mydata$position+mydata$educ_level+mydata$subject1+mydata$subject3+mydata$lang_e+mydata$genderF+mydata$age, na.action=na.exclude)
summary(model1N)
AIC(model1N)
summary(model2N)
AIC(model2N)
summary(model3N)
AIC(model3N)
summary(model4N)
AIC(model4N)

#both types of statement items, i.e. DemGov4#
model1C<-lm(mydata$DemGov4~mydata$media+mydata$Alt+mydata$polit+mydata$position+mydata$educ_level+
mydata$subject1+mydata$subject3+mydata$lang_e+mydata$genderF+mydata$age, na.action=na.exclude)
model2C<-lm(mydata$DemGov4~mydata$abroad+mydata$Alt+mydata$polit+mydata$position+mydata$educ_level+
mydata$subject1+mydata$subject3+mydata$lang_e+mydata$genderF+mydata$age, na.action=na.exclude)
model3C<-lm(mydata$DemGov4~mydata$netw+mydata$Alt+mydata$polit+mydata$position+mydata$educ_level+
mydata$subject1+mydata$subject3+mydata$lang_e+mydata$genderF+mydata$age, na.action=na.exclude)
model4C<-lm(mydata$DemGov4~mydata$media+mydata$abroad+mydata$netw+mydata$Alt+mydata$polit+
mydata$position+mydata$educ_level+mydata$subject1+mydata$subject3+mydata$lang_e+mydata$genderF+mydata$age, na.action=na.exclude)
summary(model1C)
AIC(model1C)
summary(model2C)
AIC(model2C)
summary(model3C)
AIC(model3C)
summary(model4C)
AIC(model4C)


#Figure 2. Interaction effects#
#cross-product terms transnational influences
mydata$NetAbr<-mydata$netw*mydata$abroad
mydata$NetMed<-mydata$netw*mydata$media
mydata$AbrMed<-mydata$abroad*mydata$media
model6C<-lm(mydata$DemGov4~mydata$netw +mydata$NetAbr +mydata$abroad +mydata$media+mydata$Alt+mydata$polit+mydata$position+mydata$educ_level+mydata$subject1+mydata$subject3+mydata$lang_e+mydata$genderF+mydata$age, na.action=na.exclude)
model7C<-lm(mydata$DemGov4~mydata$netw +mydata$NetMed +mydata$abroad +mydata$media+mydata$Alt+mydata$polit+mydata$position+mydata$educ_level+mydata$subject1+mydata$subject3+mydata$lang_e+mydata$genderF+mydata$age, na.action=na.exclude)
model8C<-lm(mydata$DemGov4~mydata$netw +mydata$AbrMed +mydata$abroad +mydata$media+mydata$Alt+mydata$polit+mydata$position+mydata$educ_level+mydata$subject1+mydata$subject3+mydata$lang_e+mydata$genderF+mydata$age, na.action=na.exclude)
summary(model6C)
AIC(model6C)
summary(model7C)
AIC(model7C)
summary(model8C)
AIC(model8C)

mydataPOL <- subset(mydata, polit==1)
mydataPOL$NetAbr<-mydataPOL$netw*mydataPOL$abroad
mydataPOL$NetMed<-mydataPOL$netw*mydataPOL$media
mydataPOL$AbrMed<-mydataPOL$abroad*mydataPOL$media
modelP6C<-lm(mydataPOL$DemGov4~mydataPOL$netw +mydataPOL$NetAbr +mydataPOL$abroad+mydataPOL$media+mydataPOL$Alt+mydataPOL$position+mydataPOL$educ_level+mydataPOL$subject1+mydataPOL$subject3+mydataPOL$lang_e+mydataPOL$genderF+mydataPOL$age, na.action=na.exclude)
modelP7C<-lm(mydataPOL$DemGov4~mydataPOL$netw +mydataPOL$NetMed +mydataPOL$abroad +mydataPOL$media+mydataPOL$Alt+mydataPOL$position+mydataPOL$educ_level+mydataPOL$subject1+mydataPOL$subject3+mydataPOL$lang_e+mydataPOL$genderF+mydataPOL$age, na.action=na.exclude)
modelP8C<-lm(mydataPOL$DemGov4~mydataPOL$netw +mydataPOL$AbrMed +mydataPOL$abroad +mydataPOL$media+mydataPOL$Alt+mydataPOL$position+mydataPOL$educ_level+mydataPOL$subject1+mydataPOL$subject3+mydataPOL$lang_e+mydataPOL$genderF+mydataPOL$age, na.action=na.exclude)
summary(modelP6C)
AIC(modelP6C)
summary(modelP7C)
AIC(modelP7C)
summary(modelP8C)
AIC(modelP8C)


mydata$NetAbr<-mydata$netw*mydata$abroad
mydata$NetMed<-mydata$netw*mydata$media
mydata$AbrMed<-mydata$abroad*mydata$media
model6C<-lm(mydata$DemGov1~mydata$netw +mydata$NetAbr +mydata$abroad +mydata$media+mydata$Alt+mydata$polit+mydata$position+mydata$educ_level+mydata$subject1+mydata$subject3+mydata$lang_e+mydata$genderF+mydata$age, na.action=na.exclude)
model7C<-lm(mydata$DemGov1~mydata$netw +mydata$NetMed +mydata$abroad +mydata$media+mydata$Alt+mydata$polit+mydata$position+mydata$educ_level+mydata$subject1+mydata$subject3+mydata$lang_e+mydata$genderF+mydata$age, na.action=na.exclude)
model8C<-lm(mydata$DemGov1~mydata$netw +mydata$AbrMed +mydata$abroad +mydata$media+mydata$Alt+mydata$polit+mydata$position+mydata$educ_level+mydata$subject1+mydata$subject3+mydata$lang_e+mydata$genderF+mydata$age, na.action=na.exclude)
summary(model6C)
AIC(model6C)
summary(model7C)
AIC(model7C)
summary(model8C)
AIC(model8C)

mydataPOL <- subset(mydata, polit==1)
mydataPOL$NetAbr<-mydataPOL$netw*mydataPOL$abroad
mydataPOL$NetMed<-mydataPOL$netw*mydataPOL$media
mydataPOL$AbrMed<-mydataPOL$abroad*mydataPOL$media
modelP6C<-lm(mydataPOL$DemGov1~mydataPOL$netw +mydataPOL$NetAbr +mydataPOL$abroad+mydataPOL$media+mydataPOL$Alt+mydataPOL$position+mydataPOL$educ_level+mydataPOL$subject1+mydataPOL$subject3+mydataPOL$lang_e+mydataPOL$genderF+mydataPOL$age, na.action=na.exclude)
modelP7C<-lm(mydataPOL$DemGov1~mydataPOL$netw +mydataPOL$NetMed +mydataPOL$abroad +mydataPOL$media+mydataPOL$Alt+mydataPOL$position+mydataPOL$educ_level+mydataPOL$subject1+mydataPOL$subject3+mydataPOL$lang_e+mydataPOL$genderF+mydataPOL$age, na.action=na.exclude)
modelP8C<-lm(mydataPOL$DemGov1~mydataPOL$netw +mydataPOL$AbrMed +mydataPOL$abroad +mydataPOL$media+mydataPOL$Alt+mydataPOL$position+mydataPOL$educ_level+mydataPOL$subject1+mydataPOL$subject3+mydataPOL$lang_e+mydataPOL$genderF+mydataPOL$age, na.action=na.exclude)
summary(modelP6C)
AIC(modelP6C)
summary(modelP7C)
AIC(modelP7C)
summary(modelP8C)
AIC(modelP8C)

mydata$NetAbr<-mydata$netw*mydata$abroad
mydata$NetMed<-mydata$netw*mydata$media
mydata$AbrMed<-mydata$abroad*mydata$media
model6C<-lm(mydata$DemGov3~mydata$netw +mydata$NetAbr +mydata$abroad +mydata$media+mydata$Alt+mydata$polit+mydata$position+mydata$educ_level+mydata$subject1+mydata$subject3+mydata$lang_e+mydata$genderF+mydata$age, na.action=na.exclude)
model7C<-lm(mydata$DemGov3~mydata$netw +mydata$NetMed +mydata$abroad +mydata$media+mydata$Alt+mydata$polit+mydata$position+mydata$educ_level+mydata$subject1+mydata$subject3+mydata$lang_e+mydata$genderF+mydata$age, na.action=na.exclude)
model8C<-lm(mydata$DemGov3~mydata$netw +mydata$AbrMed +mydata$abroad +mydata$media+mydata$Alt+mydata$polit+mydata$position+mydata$educ_level+mydata$subject1+mydata$subject3+mydata$lang_e+mydata$genderF+mydata$age, na.action=na.exclude)
summary(model6C)
AIC(model6C)
summary(model7C)
AIC(model7C)
summary(model8C)
AIC(model8C)

mydataPOL <- subset(mydata, polit==1)
mydataPOL$NetAbr<-mydataPOL$netw*mydataPOL$abroad
mydataPOL$NetMed<-mydataPOL$netw*mydataPOL$media
mydataPOL$AbrMed<-mydataPOL$abroad*mydataPOL$media
modelP6C<-lm(mydataPOL$DemGov3~mydataPOL$netw +mydataPOL$NetAbr +mydataPOL$abroad+mydataPOL$media+mydataPOL$Alt+mydataPOL$position+mydataPOL$educ_level+mydataPOL$subject1+mydataPOL$subject3+mydataPOL$lang_e+mydataPOL$genderF+mydataPOL$age, na.action=na.exclude)
modelP7C<-lm(mydataPOL$DemGov3~mydataPOL$netw +mydataPOL$NetMed +mydataPOL$abroad +mydataPOL$media+mydataPOL$Alt+mydataPOL$position+mydataPOL$educ_level+mydataPOL$subject1+mydataPOL$subject3+mydataPOL$lang_e+mydataPOL$genderF+mydataPOL$age, na.action=na.exclude)
modelP8C<-lm(mydataPOL$DemGov3~mydataPOL$netw +mydataPOL$AbrMed +mydataPOL$abroad +mydataPOL$media+mydataPOL$Alt+mydataPOL$position+mydataPOL$educ_level+mydataPOL$subject1+mydataPOL$subject3+mydataPOL$lang_e+mydataPOL$genderF+mydataPOL$age, na.action=na.exclude)
summary(modelP6C)
AIC(modelP6C)
summary(modelP7C)
AIC(modelP7C)
summary(modelP8C)
AIC(modelP8C)


#####################
#Multiple imputation#
#####################
mydata<-read.table("Data.csv", header=TRUE, sep=",")

##Identification of missing mechanism##

mydata[is.na(mydata)]<-99

mydata$missing <-with(mydata, ifelse (
abroad_EU==99|
abroad_USA==99 |
media==99|
age==99|
educ_level==99|
item1==99|
item2==99|
item3==99|
item4==99|
item5==99|
item6==99|
item7==99|
item8==99|
itemN3==99|
itemN1==99|
itemN2==99, 1, 0))

table(mydata$missing)
na<-density(mydata$missing)
plot(na,type="n", main="Missing Values")
polygon(na, col="lightgray", border="gray")

#Preparation of dataset##
data.items<-read.table("Items.csv", header=TRUE, sep=";")
names(data.items)
#Q22a, Q22d, Q23a, Q23b, Q23c are additional statement items
items <- c("id","Q22a","Q22d","Q23a","Q23b","Q23c")
data.items2 <- data.items[items]
mydata<-read.table("Data.csv", header=TRUE, sep=",")
names(mydata)
# exclude variables 
vars <- names(mydata) %in% c("X", "influences", "subject1", "subject2", "subject3") 
mydata <- mydata[!vars]
names(mydata)
mydata <- merge(mydata,data.items2,by="id")
head(mydata)

##imputation of missing values##
library(foreign)
library(RcppArmadillo)
library(Amelia)

vars_to_impute<-mydata[,c("ministry","department","netw_Twin","netw_EUMS","netw_UNDP",
"netw_USA","netw_WB","netw_Japon","media","abroad_EU","abroad_USA","age",
"genderF","position","educ_level","lang_e","educ_subject",
"item1","item2","item3","item4","item5","item6","item7","item8","itemN3","itemN1",
"itemN2","Q22a","Q22d","Q23a","Q23b","Q23c")]
a.out<-amelia(vars_to_impute, ords=c("ministry","department","netw_Twin","netw_EUMS","netw_UNDP",
"netw_USA","netw_WB","netw_Japon","media","abroad_EU","abroad_USA","age",
"genderF","position","educ_level","lang_e","educ_subject",
"item1","item2","item3","item4","item5","item6","item7","item8","itemN3","itemN1",
"itemN2","Q22a","Q22d","Q23a","Q23b","Q23c"))
table(round(a.out$imputations[[3]]$item1, digits=3))
summary(a.out)
save(a.out, file = "imputations.RData")
write.amelia(obj=a.out, file.stem="inmi", format="csv", na="NA")
#missmap(a.out)


##transformation of variables##
a.out <- transform(a.out, DemGov1=item1+item2+item3+item4+item5+item6+item7+item8)
a.out <- transform(a.out, DemGov3=itemN1+itemN2+itemN3)
a.out <- transform(a.out, DemGov4=DemGov1+DemGov3)
a.out <- transform(a.out, DemGov2=item1+item2+item3+item4+item5+item6+item7+item8+itemN1+itemN2+itemN3)
a.out <- transform(a.out, netw=ifelse(ministry==11,0,netw_Twin))
#a.out <- transform(a.out, Alt=ifelse(netw_EUMS==1|netw_UNDP==1|netw_USA==1|netw_WB==1|netw_Japon==1,1,0))
a.out <- transform(a.out, Alt=ifelse(netw_Twin==0&(netw_EUMS==1|netw_UNDP==1|netw_USA==1|netw_WB==1|netw_Japon==1),1,0))
a.out <- transform(a.out, abroad=ifelse(abroad_EU==1|abroad_USA==1,1,0))
a.out <- transform(a.out, infl=ifelse(abroad==1|media==1|netw==1,1,0))
a.out <- transform(a.out, polit=ifelse(ministry==2|ministry==7|ministry==8|ministry==9|ministry==10|ministry==11,1,0))
#subject1=law/econ; subject2 = natural sciences; subject3 = public admin
a.out <- transform(a.out, subject1=ifelse(educ_subject==1,1,0))
a.out <- transform(a.out, subject2=ifelse(educ_subject==2,1,0))
a.out <- transform(a.out, subject3=ifelse(educ_subject==2,1,0))
head(a.out$imputations[[1]])
summary(a.out)
save(a.out, file = "imputations.RData")
write.amelia(obj=a.out, file.stem="inmi", format="csv", na="NA")

plot(a.out, which.vars=35)

inmi1<-read.table("inmi1.csv", header=TRUE, sep=",")
names(inmi1)

################
#Selection bias#
################

model.media <- zelig(media~ministry, model="ls", data=a.out$imputations)
model.study <- zelig(abroad~ministry, model="ls", data=a.out$imputations)
model.media2 <- zelig(media~department, model="ls", data=a.out$imputations)
model.study2 <- zelig(abroad~department, model="ls", data=a.out$imputations)

media <- coef(summary(model.media)) 
print(media, digits = 2)
study <- coef(summary(model.study)) 
print(study, digits = 2)
media2 <- coef(summary(model.media2)) 
print(media2, digits = 2)
study2 <- coef(summary(model.study2)) 
print(study2, digits = 2)

table(inmi1$netw,inmi2$polit) 


t.test(a.out$imputations[[1]]$ministry~a.out$imputations[[1]]$media) # where y is numeric and x is a binary factor 
t.test(a.out$imputations[[5]]$ministry~a.out$imputations[[5]]$media)

t.test(a.out$imputations[[1]]$ministry~a.out$imputations[[1]]$abroad)
t.test(a.out$imputations[[5]]$ministry~a.out$imputations[[5]]$abroad)

t.test(a.out$imputations[[1]]$department~a.out$imputations[[1]]$media)
t.test(a.out$imputations[[5]]$department~a.out$imputations[[5]]$media)

t.test(a.out$imputations[[1]]$department~a.out$imputations[[1]]$abroad)
t.test(a.out$imputations[[5]]$department~a.out$imputations[[5]]$abroad) 


###########
#Causation#
###########
library("Zelig")
engl.N1<- zelig(lang_e~netw,model="ls",data=a.out$imputations)
engl.A1<- zelig(lang_e~abroad,model="ls",data=a.out$imputations)
posit.A2<- zelig(position~abroad,model="ls",data=a.out$imputations)
engl.N2<- zelig(lang_e~netw+educ_level+subject1+subject3+genderF+age,model="ls",data=a.out$imputations)
engl.A3<- zelig(lang_e~abroad+educ_level+subject1+subject3+genderF+age,model="ls",data=a.out$imputations)
posit.A4<- zelig(position~abroad+educ_level+subject1+subject3+genderF+age,model="ls",data=a.out$imputations)
summary(engl.N1)
summary(engl.A1)
summary(posit.A2)
summary(engl.N2)
summary(engl.A3)
summary(posit.A4)

engl.N <- xtabs(~lang_e+netw,data=a.out$imputations[[1]])
engl.A <- xtabs(~lang_e+abroad,data=a.out$imputations[[1]])
posit.A<- xtabs(~position+abroad,data=a.out$imputations[[1]])
ftable(engl.N)
summary(engl.N)
ftable(engl.A)
summary(engl.A)
ftable(posit.A)
summary(posit.A)

#####################
#Regression Analyses#
#####################
library(Zelig)

inmi1<-read.table("inmi1.csv", header=TRUE, sep=",")
inmi2<-read.table("inmi1.csv", header=TRUE, sep=",")
inmi3<-read.table("inmi1.csv", header=TRUE, sep=",")
inmi4<-read.table("inmi1.csv", header=TRUE, sep=",")
inmi5<-read.table("inmi1.csv", header=TRUE, sep=",")
names(inmi1)

#Independent effects, combined dataset Table OA.2#
model.C1 <- zelig(DemGov4~media+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age, 
model="ls", data=mi(inmi1,inmi2,inmi3,inmi4,inmi5))
model.C2 <- zelig(DemGov4~abroad+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age, 
model="ls", data=mi(inmi1,inmi2,inmi3,inmi4,inmi5))
model.C3 <- zelig(DemGov4~netw+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age,  
model="ls", data=mi(inmi1,inmi2,inmi3,inmi4,inmi5))
model.C4 <- zelig(DemGov4~media+abroad+netw+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age, 
model="ls", data=mi(inmi1,inmi2,inmi3,inmi4,inmi5))
model.P1 <- zelig(DemGov1~media+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age,
model="ls", data=mi(inmi1,inmi2,inmi3,inmi4,inmi5))
model.P2 <- zelig(DemGov1~abroad+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age, 
model="ls", data=mi(inmi1,inmi2,inmi3,inmi4,inmi5))
model.P3 <- zelig(DemGov1~netw+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age,  
model="ls", data=mi(inmi1,inmi2,inmi3,inmi4,inmi5))
model.P4 <- zelig(DemGov1~media+abroad+netw+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age, 
model="ls", data=mi(inmi1,inmi2,inmi3,inmi4,inmi5))
model.N1 <- zelig(DemGov3~media+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age,  
model="ls", data=mi(inmi1,inmi2,inmi3,inmi4,inmi5))
model.N2 <- zelig(DemGov3~abroad+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age,  
model="ls", data=mi(inmi1,inmi2,inmi3,inmi4,inmi5))
model.N3 <- zelig(DemGov3~netw+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age,   
model="ls",data=mi(inmi1,inmi2,inmi3,inmi4,inmi5))
model.N4 <- zelig(DemGov3~media+abroad+netw+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age,  
model="ls", data=mi(inmi1,inmi2,inmi3,inmi4,inmi5))

summary(model.C4) 
summary(model.P4) 
summary(model.N4) 

C1 <- summary(model.C4[[1]])
C2 <- summary(model.C4[[2]])
C3 <- summary(model.C4[[3]])
C4 <- summary(model.C4[[4]])
C5 <- summary(model.C4[[5]])
P1 <- summary(model.P4[[1]])
P2 <- summary(model.P4[[2]])
P3 <- summary(model.P4[[3]])
P4 <- summary(model.P4[[4]])
P5 <- summary(model.P4[[5]])
N1 <- summary(model.N4[[1]])
N2 <- summary(model.N4[[2]])
N3 <- summary(model.N4[[3]])
N4 <- summary(model.N4[[4]])
N5 <- summary(model.N4[[5]])

(C1$fstatistic[1]+C2$fstatistic[1]+C3$fstatistic[1]+C4$fstatistic[1]+C5$fstatistic[1])/5
(P1$fstatistic[1]+P2$fstatistic[1]+P3$fstatistic[1]+P4$fstatistic[1]+P5$fstatistic[1])/5
(N1$fstatistic[1]+N2$fstatistic[1]+N3$fstatistic[1]+N4$fstatistic[1]+N5$fstatistic[1])/5

(AIC(model.C4[[1]])+AIC(model.C4[[2]])+AIC(model.C4[[3]])+AIC(model.C4[[4]])+AIC(model.C4[[5]]))/5
(AIC(model.P4[[1]])+AIC(model.P4[[2]])+AIC(model.P4[[3]])+AIC(model.P4[[4]])+AIC(model.P4[[5]]))/5
(AIC(model.N4[[1]])+AIC(model.N4[[2]])+AIC(model.N4[[3]])+AIC(model.N4[[4]])+AIC(model.N4[[5]]))/5

########################
#Robust standard errors#
########################
library(sandwich)

inmi1<-read.table("inmi1.csv", header=TRUE, sep=",")
inmi2<-read.table("inmi1.csv", header=TRUE, sep=",")
inmi3<-read.table("inmi1.csv", header=TRUE, sep=",")
inmi4<-read.table("inmi1.csv", header=TRUE, sep=",")
inmi5<-read.table("inmi1.csv", header=TRUE, sep=",")
names(inmi1)

#Independent effects, combined dataset Table OA.2#
model.C1.r <- zelig(DemGov4~media+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age, 
model="ls", data=mi(inmi1,inmi2,inmi3,inmi4,inmi5), robust=TRUE, cluster="department")
model.C2.r <- zelig(DemGov4~abroad+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age, 
model="ls", data=mi(inmi1,inmi2,inmi3,inmi4,inmi5), robust=TRUE, cluster="department")
model.C3.r <- zelig(DemGov4~netw+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age,  
model="ls", data=mi(inmi1,inmi2,inmi3,inmi4,inmi5), robust=TRUE, cluster="department")
model.C4.r <- zelig(DemGov4~media+abroad+netw+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age, 
model="ls", data=mi(inmi1,inmi2,inmi3,inmi4,inmi5), robust=TRUE, cluster="department")
model.P1.r <- zelig(DemGov1~media+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age,
model="ls", data=mi(inmi1,inmi2,inmi3,inmi4,inmi5), robust=TRUE, cluster="department")
model.P2.r <- zelig(DemGov1~abroad+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age, 
model="ls", data=mi(inmi1,inmi2,inmi3,inmi4,inmi5), robust=TRUE, cluster="department")
model.P3.r <- zelig(DemGov1~netw+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age,  
model="ls", data=mi(inmi1,inmi2,inmi3,inmi4,inmi5), robust=TRUE, cluster="department")
model.P4.r <- zelig(DemGov1~media+abroad+netw+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age, 
model="ls", data=mi(inmi1,inmi2,inmi3,inmi4,inmi5), robust=TRUE, cluster="department")
model.N1.r <- zelig(DemGov3~media+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age,  
model="ls", data=mi(inmi1,inmi2,inmi3,inmi4,inmi5), robust=TRUE, cluster="department")
model.N2.r <- zelig(DemGov3~abroad+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age,  
model="ls", data=mi(inmi1,inmi2,inmi3,inmi4,inmi5), robust=TRUE, cluster="department")
model.N3.r <- zelig(DemGov3~netw+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age,   
model="ls",data=mi(inmi1,inmi2,inmi3,inmi4,inmi5), robust=TRUE, cluster="department")
model.N4.r <- zelig(DemGov3~media+abroad+netw+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age,  
model="ls", data=mi(inmi1,inmi2,inmi3,inmi4,inmi5), robust=TRUE, cluster="department")

summary(model.C4.r) 
summary(model.P4.r) 
summary(model.N4.r) 

C1 <- summary(model.C4.r[[1]])
C2 <- summary(model.C4.r[[2]])
C3 <- summary(model.C4.r[[3]])
C4 <- summary(model.C4.r[[4]])
C5 <- summary(model.C4.r[[5]])
P1 <- summary(model.P4.r[[1]])
P2 <- summary(model.P4.r[[2]])
P3 <- summary(model.P4.r[[3]])
P4 <- summary(model.P4.r[[4]])
P5 <- summary(model.P4.r[[5]])
N1 <- summary(model.N4.r[[1]])
N2 <- summary(model.N4.r[[2]])
N3 <- summary(model.N4.r[[3]])
N4 <- summary(model.N4.r[[4]])
N5 <- summary(model.N4.r[[5]])

(C1$fstatistic[1]+C2$fstatistic[1]+C3$fstatistic[1]+C4$fstatistic[1]+C5$fstatistic[1])/5
(P1$fstatistic[1]+P2$fstatistic[1]+P3$fstatistic[1]+P4$fstatistic[1]+P5$fstatistic[1])/5
(N1$fstatistic[1]+N2$fstatistic[1]+N3$fstatistic[1]+N4$fstatistic[1]+N5$fstatistic[1])/5

(C1$df[1]+C2$df[1]+C3$df[1]+C4$df[1]+C5$df[1])/5
(P1$df[1]+P2$df[1]+P3$df[1]+P4$df[1]+P5$df[1])/5
(N1$df[1]+N2$df[1]+N3$df[1]+N4$df[1]+N5$df[1])/5

(AIC(model.C4.r[[1]])+AIC(model.C4.r[[2]])+AIC(model.C4.r[[3]])+AIC(model.C4.r[[4]])+AIC(model.C4.r[[5]]))/5
(AIC(model.P4.r[[1]])+AIC(model.P4.r[[2]])+AIC(model.P4.r[[3]])+AIC(model.P4.r[[4]])+AIC(model.P4.r[[5]]))/5
(AIC(model.N4.r[[1]])+AIC(model.N4.r[[2]])+AIC(model.N4.r[[3]])+AIC(model.N4.r[[4]])+AIC(model.N4.r[[5]]))/5


#Figure 2. Regression Plot full model (no. 4)
var.names <- c("Foreign media", "Int. education", "Transgov. network", "Politicization", "Alternative program", "Directorate", "Postgraduate", "Law/Econ.", "Public Admin.", "English", "Women", "Age")
x1 <- summary(model.C4.r) 
x2 <- summary(model.P4.r) 
x3 <- summary(model.N4.r)

coef.vec.1 <- x1$coefficients[2:13,1]  #combined scale
se.vec.1   <- x1$coefficients[2:13,2] 
coef.vec.2<- x2$coefficients[2:13,1]  #positive
se.vec.2 <- x2$coefficients[2:13,2]
coef.vec.3<- x3$coefficients[2:13,1]  #negative
se.vec.3<-x3$coefficients[2:13,2]

pdf("Figure2aR.pdf", height = 8, width = 7)
y.axis <- length(var.names):1
adjust <- 0.1 #create object to adjust points and lines up and down to distinguish between models
#par(mfrow=c(1,2), oma=c(1,1,1,1))
par(mar=c(3,9,1,3))

plot(coef.vec.1, y.axis+adjust, type = "p", axes = F, 
xlab = list(label="Standard. regression coefficient", cex=0.9), 
ylab = list(label="", cex=0.9), pch = 19, cex = 1.2, 
xlim = c(min((coef.vec.1-qnorm(.95)*se.vec.1 -.1)),
max((coef.vec.1+qnorm(.95)*se.vec.1 -.1))),  
ylim = c(min(y.axis), max(y.axis)), main = "")
axis(1,at = seq(-3.0,5.0, by = 1), label = seq(-3.0,5.0, by = 1), cex.axis = 1.1)
axis(2, at = y.axis, label = var.names, las = 1, tick = T, cex.axis =1.1)
abline(h = y.axis, lty = 2, lwd = .5, col = "light grey")
segments(coef.vec.1-qnorm(.95)*se.vec.1, y.axis+adjust, coef.vec.1+qnorm(.95)*se.vec.1, y.axis+adjust, lwd =  1.3)
abline(v=0, lty = 2)
dev.off()

pdf("Figure2bR.pdf", height = 8, width = 7)
par(mar=c(3,9,1,3))
plot(coef.vec.2, y.axis+adjust, type = "p", axes = F, xlab = list(label="Standard. regression coefficient", cex=0.9), ylab = list(label="", cex=0.9), pch = 21, cex = 1.2, 
xlim = c(min((coef.vec.2-qnorm(.95)*se.vec.2 -.1), (coef.vec.3-qnorm(.95)*se.vec.3 -.1)), #set xlims at mins and maximums (from models) of confidence intervals, plus .1 to leave room at ends of plots
    max((coef.vec.2+qnorm(.95)*se.vec.2 -.1), (coef.vec.3+qnorm(.95)*se.vec.3 -.1))),  
ylim = c(min(y.axis), max(y.axis)), main = "")
axis(1,at = seq(-3.0,5.0, by = 1), label = seq(-3.0,5.0, by = 1), cex.axis = 1.1)
axis(2, at = y.axis,label = var.names, las = 1, tick = T, cex.axis =1.1)
abline(h = y.axis, lty = 2, lwd = .5, col = "light grey")
segments(coef.vec.2-qnorm(.95)*se.vec.2, y.axis+adjust, coef.vec.2+qnorm(.95)*se.vec.2, y.axis+adjust, lwd =  1.3)
abline(v=0, lty = 2)
#add models
segments(coef.vec.3-qnorm(.95)*se.vec.3, y.axis-adjust, coef.vec.3+qnorm(.95)*se.vec.3, y.axis-adjust, lwd =  1.3)
points(coef.vec.3, y.axis-adjust, pch = 17, cex = 1.2, bg = "black" )
#legend(xpd=TRUE, -11,-1,c("Combined scale", "Positively-oriented scale", "Negatively-oriented scale"), pch = c(19,21,17),lty = c(1,3))
legend("bottomleft",border=TRUE, c("posit. items", "negat. items"), pch = c(21,17), cex=0.9)
dev.off()


###################
#Multi-level model#
###################

library(lattice)
inmi1<-read.table("inmi1.csv", header=TRUE, sep=",")
inmi2<-read.table("inmi2.csv", header=TRUE, sep=",")
inmi3<-read.table("inmi3.csv", header=TRUE, sep=",")
inmi4<-read.table("inmi4.csv", header=TRUE, sep=",")
inmi5<-read.table("inmi5.csv", header=TRUE, sep=",")

site.1 <- with(inmi1, reorder(department, DemGov2, FUN=mean))
stripplot(site.1 ~ DemGov2, data=inmi1, scales = list(tck=0.5))
site.2 <- with(inmi2, reorder(department, DemGov2, FUN=mean))
stripplot(site.2 ~ DemGov2, data=inmi1, scales = list(tck=0.5))
site.3 <- with(inmi3, reorder(department, DemGov2, FUN=mean))
stripplot(site.3 ~ DemGov2, data=inmi1, scales = list(tck=0.5))
site.4 <- with(inmi4, reorder(department, DemGov2, FUN=mean))
stripplot(site.4 ~ DemGov2, data=inmi1, scales = list(tck=0.5))
site.5 <- with(inmi5, reorder(department, DemGov2, FUN=mean))
stripplot(site.5 ~ DemGov2, data=inmi1, scales = list(tck=0.5))

library(Zelig)
library(lme4) #https://stat.ethz.ch/pipermail/r-help/2006-May/094765.html
library(languageR)

#null-model
model.0 <- zelig(DemGov4~1, model="ls", data=a.out$imputations)
m.0 <- coef(summary(model.0)) 
print(m.0, digits = 2)

ML.model.0 <-list(5)  # a list to store each model
for (i in 1:5) {
file.name <- paste("inmi", 5 ,".csv",sep="")
data.to.use <- read.csv(file.name)
ML.model.0[[5]] <- lmer(DemGov4 ~ 1 + (1 | department),
data = data.to.use)}
m.0 <- ML.model.0[[5]] 
pvals.fnc(m.0)
ranef(m.0)

#random intercept, fixed predictor in individual level, REML = residual maximum likelihood
ML.model.C4 <-list(5)  # a list to store each model
for (i in 1:5) {
file.name <- paste("inmi", 5 ,".csv",sep="")
data.to.use <- read.csv(file.name)
ML.model.C4[[5]] <- lmer(DemGov4 ~ media+abroad+netw+Alt+position+polit+educ_level+subject1+subject3+lang_e+genderF+age +
(1| department), data = data.to.use)}
m.C4 <- ML.model.C4[[5]] 
pvals.fnc(m.C4)
ranef(m.C4)

ML.model.C2 <-list(5)  # a list to store each model
for (i in 1:5) {
file.name <- paste("inmi", 5 ,".csv",sep="")
data.to.use <- read.csv(file.name)
ML.model.C2[[5]] <- lmer(DemGov2 ~ media+abroad+netw+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age +
(1| department), data = data.to.use)}
m.C2 <- ML.model.C2[[5]] 
pvals.fnc(m.C2)
ranef(m.C2)

ML.model.P <-list(5) 
for (i in 1:5) {
file.name <- paste("inmi", 5 ,".csv",sep="")
data.to.use <- read.csv(file.name)
ML.model.P[[5]] <- lmer(DemGov1 ~ media+abroad+netw+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age + 
(1 | department), data = data.to.use)}
m.P <- ML.model.P[[5]] 
pvals.fnc(m.P)
ranef(m.P)

ML.model.N <-list(5)  
for (i in 1:5) {
file.name <- paste("inmi", 5 ,".csv",sep="")
data.to.use <- read.csv(file.name)
ML.model.N[[5]] <- lmer(DemGov3 ~ media+abroad+netw+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age +
(1 | department), data = data.to.use)}
m.N <- ML.model.N[[5]] 
pvals.fnc(m.N)
ranef(m.N)

#####################
#INTERACTION EFFECTS#
#####################

#Regression on imputed data
library(Zelig)
#a.out$imputations muss im memory sein! here with robust standard errors

#two-way interaction: transnational influences
model.C5 <- zelig(DemGov4~media+(netw*abroad)+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age, 
model="ls", data=a.out$imputations, robust=TRUE, cluster="department")
model.C6 <- zelig(DemGov4~abroad+(netw*media)+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age, 
model="ls", data=a.out$imputations, robust=TRUE, cluster="department")
model.P5 <- zelig(DemGov1~media+(netw*abroad)+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age, 
model="ls", data=a.out$imputations, robust=TRUE, cluster="department")
model.P6 <- zelig(DemGov1~abroad+(netw*media)+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age, 
model="ls", data=a.out$imputations, robust=TRUE, cluster="department")
model.N5 <- zelig(DemGov3~media+(netw*abroad)+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age, 
model="ls", data=a.out$imputations, robust=TRUE, cluster="department")
model.N6 <- zelig(DemGov3~abroad+(netw*media)+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age, 
model="ls", data=a.out$imputations, robust=TRUE, cluster="department")
summary(model.C5)
summary(model.C6)
summary(model.P5)
summary(model.P6)
summary(model.N5)
summary(model.N6)

C5.1 <- summary(model.C5[[1]])
C5.2 <- summary(model.C5[[2]])
C5.3 <- summary(model.C5[[3]])
C5.4 <- summary(model.C5[[4]])
C5.5 <- summary(model.C5[[5]])
P5.1 <- summary(model.P5[[1]])
P5.2 <- summary(model.P5[[2]])
P5.3 <- summary(model.P5[[3]])
P5.4 <- summary(model.P5[[4]])
P5.5 <- summary(model.P5[[5]])
N5.1 <- summary(model.N5[[1]])
N5.2 <- summary(model.N5[[2]])
N5.3 <- summary(model.N5[[3]])
N5.4 <- summary(model.N5[[4]])
N5.5 <- summary(model.N5[[5]])

C6.1 <- summary(model.C6[[1]])
C6.2 <- summary(model.C6[[2]])
C6.3 <- summary(model.C6[[3]])
C6.4 <- summary(model.C6[[4]])
C6.5 <- summary(model.C6[[5]])
P6.1 <- summary(model.P6[[1]])
P6.2 <- summary(model.P6[[2]])
P6.3 <- summary(model.P6[[3]])
P6.4 <- summary(model.P6[[4]])
P6.5 <- summary(model.P6[[5]])
N6.1 <- summary(model.N6[[1]])
N6.2 <- summary(model.N6[[2]])
N6.3 <- summary(model.N6[[3]])
N6.4 <- summary(model.N6[[4]])
N6.5 <- summary(model.N6[[5]])

(C5.1$fstatistic[1]+C5.2$fstatistic[1]+C5.3$fstatistic[1]+C5.4$fstatistic[1]+C5.5$fstatistic[1])/5
(C6.1$fstatistic[1]+C6.2$fstatistic[1]+C6.3$fstatistic[1]+C6.4$fstatistic[1]+C6.5$fstatistic[1])/5
(P5.1$fstatistic[1]+P5.2$fstatistic[1]+P5.3$fstatistic[1]+P5.4$fstatistic[1]+P5.5$fstatistic[1])/5
(P6.1$fstatistic[1]+P6.2$fstatistic[1]+P6.3$fstatistic[1]+P6.4$fstatistic[1]+P6.5$fstatistic[1])/5
(N5.1$fstatistic[1]+N5.2$fstatistic[1]+N5.3$fstatistic[1]+N5.4$fstatistic[1]+N5.5$fstatistic[1])/5
(N6.1$fstatistic[1]+N6.2$fstatistic[1]+N6.3$fstatistic[1]+N6.4$fstatistic[1]+N6.5$fstatistic[1])/5

(C5.1$df[1]+C5.2$df[1]+C5.3$df[1]+C5.4$df[1]+C5.5$df[1])/5
(C6.1$df[1]+C6.2$df[1]+C6.3$df[1]+C6.4$df[1]+C6.5$df[1])/5
(P5.1$df[1]+P5.2$df[1]+P5.3$df[1]+P5.4$df[1]+P5.5$df[1])/5
(P6.1$df[1]+P6.2$df[1]+P6.3$df[1]+P6.4$df[1]+P6.5$df[1])/5
(N5.1$df[1]+N5.2$df[1]+N5.3$df[1]+N5.4$df[1]+N5.5$df[1])/5
(N6.1$df[1]+N6.2$df[1]+N6.3$df[1]+N6.4$df[1]+N6.5$df[1])/5

(AIC(model.C5[[1]])+AIC(model.C5[[2]])+AIC(model.C5[[3]])+AIC(model.C5[[4]])+AIC(model.C5[[5]]))/5
(AIC(model.C6[[1]])+AIC(model.C6[[2]])+AIC(model.C6[[3]])+AIC(model.C6[[4]])+AIC(model.C6[[5]]))/5
(AIC(model.P5[[1]])+AIC(model.P5[[2]])+AIC(model.P5[[3]])+AIC(model.P5[[4]])+AIC(model.P5[[5]]))/5
(AIC(model.P6[[1]])+AIC(model.P6[[2]])+AIC(model.P6[[3]])+AIC(model.P6[[4]])+AIC(model.P6[[5]]))/5
(AIC(model.N5[[1]])+AIC(model.N5[[2]])+AIC(model.N5[[3]])+AIC(model.N5[[4]])+AIC(model.N5[[5]]))/5
(AIC(model.N6[[1]])+AIC(model.N6[[2]])+AIC(model.N6[[3]])+AIC(model.N6[[4]])+AIC(model.N6[[5]]))/5


#three-way interaction: transta. influences and politicization
model.C7a <- zelig(DemGov4~media+(netw*abroad*polit)+Alt+position+educ_level+subject1+subject3+lang_e+genderF+age, 
model="ls", data=a.out$imputations)
model.C7b <- zelig(DemGov4~abroad+(netw*media*polit)+Alt+position+educ_level+subject1+subject3+lang_e+genderF+age, 
model="ls", data=a.out$imputations)
summary(model.C7a)
summary(model.C7b)

#two-way interaction on separate data sets: When we find complex interactions, we often break down the data into smaller data sets, for ease of interpretation.
df1 <- do.call(rbind,lapply(a.out$imputations,function(x) x[x$polit == 1,]))
df2 <- do.call(rbind,lapply(a.out$imputations,function(x) x[x$polit != 1,]))
#netw*abroad
model.C6p <- zelig(DemGov4~media+(netw*abroad)+Alt+position+educ_level+subject1+subject3+lang_e+genderF+age, 
model="ls", data=df1)
model.C6np <- zelig(DemGov4~media+(netw*abroad)+Alt+position+educ_level+subject1+subject3+lang_e+genderF+age, 
model="ls", data=df2)
summary(model.C6p)
summary(model.C6np)

#DUMMIES
a.out <- transform(a.out, D1=ifelse(netw==0&abroad==0,1,0))
a.out <- transform(a.out, D2=ifelse(netw==1&abroad==0,1,0))
a.out <- transform(a.out, D3=ifelse(netw==0&abroad==1,1,0))
a.out <- transform(a.out, D4=ifelse(netw==1&abroad==1,1,0))

model.DC <- zelig(DemGov4~D2+D3+D4+Alt+position+educ_level+subject1+subject3+lang_e+genderF+age, 
model="ls", data=a.out$imputations, robust=TRUE)
model.DP <- zelig(DemGov1~D2+D3+D4+Alt+position+educ_level+subject1+subject3+lang_e+genderF+age, 
model="ls", data=a.out$imputations, robust=TRUE)
model.DN <- zelig(DemGov3~D2+D3+D4+Alt+position+educ_level+subject1+subject3+lang_e+genderF+age, 
model="ls", data=a.out$imputations, robust=TRUE)
model.DN2 <- zelig(DemGov3~D2+D3+D4, 
model="ls", data=a.out$imputations, robust=TRUE)
summary(model.DC)
summary(model.DP)
summary(model.DN)

library(car)
linearHypothesis(model.DC, "D4 = D3", test="Chisq")
linearHypothesis(model.DP, "D4 = D3", test="Chisq")
linearHypothesis(model.DN, "D4 = D3", test="Chisq")


###################
#ROBUSTNESS CHECKS#
###################

#ANOVA comparing multiply-imputed and complete-case datasets (Appendix 2)
#complete-case analyses
mydata<-read.table("Data.csv", header=TRUE, sep=",")
mydata$DemGov1<-(mydata$item1 + mydata$item2 + mydata$item3 + mydata$item4 + mydata$item5 + mydata$item6 + mydata$item7 + mydata$item8)
mydata$DemGov2<-(mydata$item1 + mydata$item2 + mydata$item3 + mydata$item5 + mydata$item6 + mydata$item7 + mydata$item8)
mydata$DemGov3<-(mydata$itemN1 + mydata$itemN2 + mydata$itemN3)
mydata$DemGov4<-(mydata$DemGov1 + mydata$DemGov3)
mydata$netw<-ifelse(mydata$netw_Twin==1,1,0)
mydata$netw <-ifelse(mydata$ministry==11,0,mydata$netw)
mydata$Alt<-ifelse(mydata$netw_Twin==0 & (mydata$netw_EUMS==1 | mydata$netw_UNDP==1 |
mydata$netw_USA==1 | mydata$netw_WB==1 | mydata$netw_Japon==1), 1, 0)
mydata$abroad <- ifelse(mydata$abroad_EU == 1 | mydata$abroad_USA == 1,1,0)
mydata$infl<-ifelse(mydata$abroad==1|mydata$netw==1|mydata$media==1,1,0)
mydata$polit <- ifelse(mydata$ministry==2 | mydata$ministry==7 | mydata$ministry==8 |
mydata$ministry==9 | mydata$ministry==10 | mydata$ministry==11,1,0)

model4P<-lm(mydata$DemGov1~mydata$media+mydata$abroad+mydata$netw+mydata$Alt+mydata$polit+
mydata$position+mydata$educ_level+mydata$subject1+mydata$subject3+mydata$lang_e+mydata$genderF+mydata$age, na.action=na.exclude)
model4N<-lm(mydata$DemGov3~mydata$media+mydata$abroad+mydata$netw+mydata$Alt+mydata$polit+
mydata$position+mydata$educ_level+mydata$subject1+mydata$subject3+mydata$lang_e+mydata$genderF+mydata$age, na.action=na.exclude)
model4C<-lm(mydata$DemGov4~mydata$media+mydata$abroad+mydata$netw+mydata$Alt+mydata$polit+
mydata$position+mydata$educ_level+mydata$subject1+mydata$subject3+mydata$lang_e+mydata$genderF+mydata$age, na.action=na.exclude)
summary(model4C)
summary(model4N)
summary(model4P)

anova(model4C, model.C4.r) #complete, imputed
anova(model4P, model.P4.r)
anova(model4N, model.N4.r)


#ROBUSTNESS CHECK 1: dropping item 4 due to low correlation with scale overall#
#Regression using scale without item 4
mydata2$DemGov4.x<-(mydata2$item1+mydata2$item2+mydata2$item3+mydata2$item5+mydata2$item6+mydata2$item7+mydata2$item8++mydata2$itemN1+mydata2$itemN2+mydata2$itemN3)
mydata2$DemGov1.x<-(mydata2$item1+mydata2$item2+mydata2$item3+mydata2$item5+mydata2$item6+mydata2$item7+mydata2$item8)
model.C4.x <- lm(DemGov4.x~media+abroad+netw+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age, data=mydata2)
model.P4.x <- lm(DemGov1.x~media+abroad+netw+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age, data=mydata2)
summary(model.C4.x)
summary(model.P4.x)

anova(model.C4.x, model.C4)
anova(model.P4.x, model.P4)



#ROBUSTNESS CHECK 2: matching data on all three treatment variables#
#do matching on each imputed data set, and then combine it in zelig() using mi() function

library(MatchIt)
library(Matching)
library(lattice)
library(gridExtra)

a.out.1<-read.table("inmi1.csv", header=TRUE, sep=",")
a.out.2<-read.table("inmi2.csv", header=TRUE, sep=",")
a.out.3<-read.table("inmi3.csv", header=TRUE, sep=",")
a.out.4<-read.table("inmi4.csv", header=TRUE, sep=",")
a.out.5<-read.table("inmi5.csv", header=TRUE, sep=",")
names(a.out.1)

library(reshape)
a.out.1 <- rename(a.out.1, c(X="id"))
a.out.2 <- rename(a.out.2, c(X="id"))
a.out.3 <- rename(a.out.3, c(X="id"))
a.out.4 <- rename(a.out.4, c(X="id"))
a.out.5 <- rename(a.out.5, c(X="id"))

a.out.1 <- a.out.1[,c("id","DemGov4","DemGov1","DemGov3","netw","abroad","media","age","genderF","lang_e","educ_level", "Alt", "position", "polit", "subject1", "subject2", "subject3")]
a.out.2 <- a.out.2[,c("id","DemGov4","DemGov1","DemGov3","netw","abroad","media","age","genderF","lang_e","educ_level", "Alt", "position", "polit", "subject1", "subject2", "subject3")]
a.out.3 <- a.out.3[,c("id","DemGov4","DemGov1","DemGov3","DemGov3","netw","abroad","media","age","genderF","lang_e","educ_level", "Alt", "position", "polit", "subject1", "subject2", "subject3")]
a.out.4 <- a.out.4[,c("id","DemGov4","DemGov1","DemGov3","netw","abroad","media","age","genderF","lang_e","educ_level", "Alt", "position", "polit", "subject1", "subject2", "subject3")]
a.out.5 <- a.out.5[,c("id","DemGov4","DemGov1","DemGov3","netw","abroad","media","age","genderF","lang_e","educ_level", "Alt", "position", "polit", "subject1", "subject2", "subject3")]

#ABROAD#
#Matching of treatment and control groups#
mydata.A1<-matchit(abroad~media+age+genderF+lang_e+educ_level, discard="both", method="genetic", data=a.out.1, pop.size=150)
mydata.A2<-matchit(abroad~media+age+genderF+lang_e+educ_level, discard="both", method="genetic", data=a.out.2, pop.size=150)
mydata.A3<-matchit(abroad~media+age+genderF+lang_e+educ_level, discard="both", method="genetic", data=a.out.3, pop.size=150)
mydata.A4<-matchit(abroad~media+age+genderF+lang_e+educ_level, discard="both", method="genetic", data=a.out.4, pop.size=150)
mydata.A5<-matchit(abroad~media+age+genderF+lang_e+educ_level, discard="both", method="genetic", data=a.out.5, pop.size=150)
summary(mydata.A1)
summary(mydata.A2)
summary(mydata.A3)
summary(mydata.A4)
summary(mydata.A5)
mydata.A1<-match.data(mydata.A1)
mydata.A2<-match.data(mydata.A2)
mydata.A3<-match.data(mydata.A3)
mydata.A4<-match.data(mydata.A4)
mydata.A5<-match.data(mydata.A5)

mydata<-read.table("Data.csv", header=TRUE, sep=",")
data.dep <- mydata[,c("id","department")]
mydata.A1 <- merge(data.dep,mydata.A1,by="id")
mydata.A2 <- merge(data.dep,mydata.A2,by="id")
mydata.A3 <- merge(data.dep,mydata.A3,by="id")
mydata.A4 <- merge(data.dep,mydata.A4,by="id")
mydata.A5 <- merge(data.dep,mydata.A5,by="id")

write.csv(mydata.A1, "DataA1.csv", na="NA")
write.csv(mydata.A2, "DataA2.csv", na="NA")
write.csv(mydata.A3, "DataA3.csv", na="NA")
write.csv(mydata.A4, "DataA4.csv", na="NA")
write.csv(mydata.A5, "DataA5.csv", na="NA")


#Checking of Covariate Balance with merged data set#
A1_balance_1<-MatchBalance(abroad~media+age+genderF+lang_e+educ_level, data=mydata.A1, ks=TRUE, nboot=10000, digits=3, print.level=1)
A1_balance_0<-MatchBalance(abroad~media+age+genderF+lang_e+educ_level, data=a.out.1, ks=TRUE, nboot=10000, digits=3, print.level=1)
A1_pre.treated<-subset(a.out.1, subset=(abroad=="1"))
A1_pre.control<-subset(a.out.1, subset=(abroad=="0"))
A1_post.treated<-subset(mydata.A1, subset=(abroad=="1"))
A1_post.control<-subset(mydata.A1, subset=(abroad=="0"))
sapply(A1_pre.treated, sd)
sapply(A1_pre.control, sd)
sapply(A1_post.treated, sd)
sapply(A1_post.control, sd)

#[... for all the other datasets!]


#MEDIA#
#Matching of treatment and control groups#
mydata.B1<-matchit(media~abroad+age+genderF+lang_e+educ_level, discard="both", method="genetic", data=a.out.1, pop.size=150)
mydata.B2<-matchit(media~abroad+age+genderF+lang_e+educ_level, discard="both", method="genetic", data=a.out.2, pop.size=150)
mydata.B3<-matchit(media~abroad+age+genderF+lang_e+educ_level, discard="both", method="genetic", data=a.out.3, pop.size=150)
mydata.B4<-matchit(media~abroad+age+genderF+lang_e+educ_level, discard="both", method="genetic", data=a.out.4, pop.size=150)
mydata.B5<-matchit(media~abroad+age+genderF+lang_e+educ_level, discard="both", method="genetic", data=a.out.5, pop.size=150)
summary(mydata.B1)
summary(mydata.B2)
summary(mydata.B3)
summary(mydata.B4)
summary(mydata.B5)
mydata.B1<-match.data(mydata.B1)
mydata.B2<-match.data(mydata.B2)
mydata.B3<-match.data(mydata.B3)
mydata.B4<-match.data(mydata.B4)
mydata.B5<-match.data(mydata.B5)

mydata.B1 <- merge(data.dep,mydata.B1,by="id")
mydata.B2 <- merge(data.dep,mydata.B2,by="id")
mydata.B3 <- merge(data.dep,mydata.B3,by="id")
mydata.B4 <- merge(data.dep,mydata.B4,by="id")
mydata.B5 <- merge(data.dep,mydata.B5,by="id")

write.csv(mydata.B1, "DataB1.csv", na="NA")
write.csv(mydata.B2, "DataB2.csv", na="NA")
write.csv(mydata.B3, "DataB3.csv", na="NA")
write.csv(mydata.B4, "DataB4.csv", na="NA")
write.csv(mydata.B5, "DataB5.csv", na="NA")

#Checking of Covariate Balance with merged data set#
#[... see above!]


#Regression
library(Zelig)
Amatch1<-read.table("DataA1.csv", header=TRUE, sep=",")
Amatch2<-read.table("DataA2.csv", header=TRUE, sep=",")
Amatch3<-read.table("DataA3.csv", header=TRUE, sep=",")
Amatch4<-read.table("DataA4.csv", header=TRUE, sep=",")
Amatch5<-read.table("DataA5.csv", header=TRUE, sep=",")
names(Amatch1)

library(reshape)
Amatch1 <- rename(Amatch1, c(department.x="department"))
Amatch2 <- rename(Amatch2, c(department.x="department"))
Amatch3 <- rename(Amatch3, c(department.x="department"))
Amatch4 <- rename(Amatch4, c(department.x="department"))
Amatch5 <- rename(Amatch5, c(department.x="department"))
Amatch1$department.y <- NULL
Amatch2$department.y <- NULL
Amatch3$department.y <- NULL
Amatch4$department.y <- NULL
Amatch5$department.y <- NULL

Bmatch1<-read.table("DataB1.csv", header=TRUE, sep=",")
Bmatch2<-read.table("DataB2.csv", header=TRUE, sep=",")
Bmatch3<-read.table("DataB3.csv", header=TRUE, sep=",")
Bmatch4<-read.table("DataB4.csv", header=TRUE, sep=",")
Bmatch5<-read.table("DataB5.csv", header=TRUE, sep=",")
names(Bmatch1)


A.model.C4 <- zelig(DemGov4~media+abroad+netw+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age,
model="ls", data=mi(Amatch1,Amatch2,Amatch3,Amatch4,Amatch5), robust=TRUE, cluster="department")
A.model.P4 <- zelig(DemGov1~media+abroad+netw+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age,
model="ls", data=mi(Amatch1,Amatch2,Amatch3,Amatch4,Amatch5), robust=TRUE, cluster="department")
A.model.N4 <- zelig(DemGov3~media+abroad+netw+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age,
model="ls", data=mi(Amatch1,Amatch2,Amatch3,Amatch4,Amatch5), robust=TRUE, cluster="department")

M.model.C4 <- zelig(DemGov4~media+abroad+netw+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age,
model="ls", data=mi(Bmatch1,Bmatch2,Bmatch3,Bmatch4,Bmatch5), robust=TRUE, cluster="department")
M.model.P4 <- zelig(DemGov1~media+abroad+netw+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age,
model="ls", data=mi(Bmatch1,Bmatch2,Bmatch3,Bmatch4,Bmatch5), robust=TRUE, cluster="department")
M.model.N4 <- zelig(DemGov3~media+abroad+netw+Alt+polit+position+educ_level+subject1+subject3+lang_e+genderF+age,
model="ls", data=mi(Bmatch1,Bmatch2,Bmatch3,Bmatch4,Bmatch5), robust=TRUE, cluster="department")

summary(A.model.C4)
summary(A.model.P4)
summary(A.model.N4)
summary(M.model.C4)
summary(M.model.P4)
summary(M.model.N4)

A1 <- summary(A.model.C4[[1]])
A2 <- summary(A.model.C4[[2]])
A3 <- summary(A.model.C4[[3]])
A4 <- summary(A.model.C4[[4]])
A5 <- summary(A.model.C4[[5]])
A6 <- summary(A.model.P4[[1]])
A7 <- summary(A.model.P4[[2]])
A8 <- summary(A.model.P4[[3]])
A9 <- summary(A.model.P4[[4]])
A10 <- summary(A.model.P4[[5]])
A11 <- summary(A.model.N4[[1]])
A12 <- summary(A.model.N4[[2]])
A13 <- summary(A.model.N4[[3]])
A14 <- summary(A.model.N4[[4]])
A15 <- summary(A.model.N4[[5]])
(A1$fstatistic[1]+A2$fstatistic[1]+A3$fstatistic[1]+A4$fstatistic[1]+A5$fstatistic[1])/5
(A6$fstatistic[1]+A7$fstatistic[1]+A8$fstatistic[1]+A9$fstatistic[1]+A10$fstatistic[1])/5
(A11$fstatistic[1]+A12$fstatistic[1]+A13$fstatistic[1]+A14$fstatistic[1]+A15$fstatistic[1])/5

(AIC(A.model.C4[[1]])+AIC(A.model.C4[[2]])+AIC(A.model.C4[[3]])+AIC(A.model.C4[[4]])+AIC(A.model.C4[[5]]))/5
(AIC(A.model.P4[[1]])+AIC(A.model.P4[[2]])+AIC(A.model.P4[[3]])+AIC(A.model.P4[[4]])+AIC(A.model.P4[[5]]))/5
(AIC(A.model.N4[[1]])+AIC(A.model.N4[[2]])+AIC(A.model.N4[[3]])+AIC(A.model.N4[[4]])+AIC(A.model.N4[[5]]))/5

M1 <- summary(M.model.C4[[1]])
M2 <- summary(M.model.C4[[2]])
M3 <- summary(M.model.C4[[3]])
M4 <- summary(M.model.C4[[4]])
M5 <- summary(M.model.C4[[5]])
M6 <- summary(M.model.P4[[1]])
M7 <- summary(M.model.P4[[2]])
M8 <- summary(M.model.P4[[3]])
M9 <- summary(M.model.P4[[4]])
M10 <- summary(M.model.P4[[5]])
M11 <- summary(M.model.N4[[1]])
M12 <- summary(M.model.N4[[2]])
M13 <- summary(M.model.N4[[3]])
M14 <- summary(M.model.N4[[4]])
M15 <- summary(M.model.N4[[5]])
(M1$fstatistic[1]+M2$fstatistic[1]+M3$fstatistic[1]+M4$fstatistic[1]+M5$fstatistic[1])/5
(M6$fstatistic[1]+M7$fstatistic[1]+M8$fstatistic[1]+M9$fstatistic[1]+M10$fstatistic[1])/5
(M11$fstatistic[1]+M12$fstatistic[1]+M13$fstatistic[1]+M14$fstatistic[1]+M15$fstatistic[1])/5

(AIC(M.model.C4[[1]])+AIC(M.model.C4[[2]])+AIC(M.model.C4[[3]])+AIC(M.model.C4[[4]])+AIC(M.model.C4[[5]]))/5
(AIC(M.model.P4[[1]])+AIC(M.model.P4[[2]])+AIC(M.model.P4[[3]])+AIC(M.model.P4[[4]])+AIC(M.model.P4[[5]]))/5
(AIC(M.model.N4[[1]])+AIC(M.model.N4[[2]])+AIC(M.model.N4[[3]])+AIC(M.model.N4[[4]])+AIC(M.model.N4[[5]]))/5

